#include "nBodySimulator.h"

// Initialize the simulation by reading in the parameters and the
// initial setup.
nBodySimulator::nBodySimulator(string parfn, string initfn)
{
   params  = new ParamHolder(parfn);
   holder = new nBodyHolder(initfn, params);
}

// Run the simulation using the specified parameters.
void nBodySimulator::run()
{
   double dt;
   OutputNumber = 0;

   t = params -> StartTime;

   double nextOutput = 0.0;

   // Advance by timesteps until the stop time is reached.
   while (t < params -> StopTime) {
      // Try the default time step.
      dt = params -> dt;

      if (nextOutput > params -> StopTime)
         nextOutput = params -> StopTime;

      // If an output occurs before that, run up to that output.
      if (nextOutput < t+dt) {
         dt = nextOutput - t;

	 advanceSim(dt);
	 t = nextOutput;
	 nextOutput += params -> OutputDt;

         // Figure out the name of the output file.
	 stringstream ss;

	 ss << "data/" << params -> OutputPrefix <<  setfill('0') << setw(4) << OutputNumber;
	 string ofName = ss.str();

	 holder -> writeToFile(ofName);

	 OutputNumber++;
         }

      // Otherwise, just run one timestep.
      else {
         advanceSim(dt);
	 t += dt;
         }

      }
}

// This is what actually update the poitions and velocities.
void nBodySimulator::advanceSim(double dt)
{
   int n = holder -> size();

   // Compute the acceleration.
   vec* accels = new vec[n];

   // Update in an eulerian fashion.
   for (int i = 0; i < n; i++)
      accels[i] = holder -> accel(i);

   for (int i = 0; i < n; i++) {
      Body* b = holder -> getBody(i);

      b -> vel[X] += accels[i][X] * dt;
      b -> vel[Y] += accels[i][Y] * dt;
      b -> vel[Z] += accels[i][Z] * dt;

      b -> pos[X] += b -> vel[X] * dt;
      b -> pos[Y] += b -> vel[Y] * dt;
      b -> pos[Z] += b -> vel[Z] * dt;
      }

   delete [] accels;
}
